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CHAPTER I 


INTRODUCTION 

"Because the measurement of the global gravity is a primary objective for both 
solid-earth dynamics and ocean dynamics and is an important secondary objective 
for continental geology, the determination of an improved gravitational field 
through space measurement should be an objective of highest priority for the 
1980's," 

— The Space Science Board, A Strategy for 
Earth Science from Space in the 1980's, the 
National Academy of Sciences, Washington, 

D.C., [1982]. 

In the early 1980's, NASA proposed a Geopotential Research Mission 
(GRM) to globally determine high precision gravitational and magnetic fields of 
Earth with advanced science techniques using full Earth coverage, polar orbiting 
satellites. These detailed fields are required for the improvement of Earth's 
mathematical models that help scientists to understand the geodynamical activities 
which are continually evolving the Earth's internal and external structure and the 
dynamics of the oceans that influence ocean circulation and global climates. In 
addition, the measurement of the geopotential to high degree and order will improve 
orbit determination of other geodectic satellites. 

Currently, three senarios exist for recovering the higher degree and order 
coefficients of Earth's gravity potential. The first senario, described as the "low- 
low" configuration, consists of two coplanar, low altitude, polar orbiting, satellites 
(Figure 1.1) that measure the Earth's gravity anomalies through range-rates derived 
from the Satellite-to-Satellite Tracking system, SST (discussed in Section 1.3) 
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[. Keating et al, 1986]. The second senario consists of a single low altitude, polar 
orbiting satellite equipped with a high precision gravity gradiometer sensor 
onboard. The third senario, described as the "high-low" configuration, consists of a 
high altitude satellite along with a low altitude, polar orbiting, satellite. This high- 
low configuration was also designed to use a Satellite- to-Satellite Tracking system 
(SST) for measuring Earth gravity anomalies [Holmes et al, 1987]. 

Each low altitude satellite in any case is "drag-free", which means that orbital 
effects due to nonconservative forces such as drag and solar radiation pressure are 
offset by means of a Disturbance Compensation System (DISCOS) that 
appropriately fires thrusters to offset the effects due to these disturbances. The 
DISCOS is discussed in more detail in Section 1.4 and Chapter 4. 

The orbital requirements of the low-low gravity mission also satisfy the 
orbital requirements for the magnetic mission. Therefore, the lead low-low satellite 
will measure the Earth's magnetic fields with the use of onboard scalar and vector 
magnetometers. 

To properly plan for this mission, a high precision (true) orbit simulation of 
the low-low GRM satellite mission was performed at the University of Texas at 
Austin using the Cray X-MP/24 supercomputer. With this orbit simulation and a 
nominal orbit simulation that was fitted to the true orbits in a least squares sense, 
the geodetic science community will be able to test gravity evaluation techniques 
that are essential to recovering a geopotential of high degree and order. In addition, 
a simulation of the DISCOS system was also performed to assist in the planning for 
the mission. Therefore, the purpose of this report is: 1) to present a simulation 
study of a proposed GRM mission by numerically integrating a high precision orbit 
for each low-low satellite using a geopotential with coefficients up to degree and 
order 360 and calculating the SST's one-way Doppler range-rate measurements, 2) 



to present a nominal orbit simulation of the aforementioned case with differences 
between the true and nominal ephemerides kept within specified limits, and 3) to 
develop a thruster control algorithm which simulates the DISCOS system so that 
fuel expenditure and thruster-on times can be studied. Since each of the 
aforementioned senarios contain a low altitude segment, the techniques used to 
analyze the low-low senario will be applicable to all three senarios. 

1.1 Description of the Low-Low Scenario for the 
Geopotential Research Mission 

As proposed by Keating et al, [1986], the GRM mission that was once 
referred to as GRAVS AT/MAGS AT, is configured using two 160 km altitude, 
polar orbiting, drag-free satellites in coincident, near circular orbits with the 
longitude of the ascending node equal to 90’, and the mean argument of periapse 
equal to 90*. In spite of perturbations due to the oblateness, the nature of the orbits 
is such that the mean argument of periapse is constant, a "frozen orbit" [White, 
1987]. The minim um mission lifetime goal was set at 6 months to achieve high 
resolution (1’ by 1*) groundtrack coverage of Earth which will result in an 
equatorial spacing of approximately 111 km . The SST onboard will use an 
integrated one-way Doppler system to calculate the line of sight range-rates between 
the satellites which will have separation distances of 100 km to 600 km. A 
Disturbance Compensation System (DISCOS) will keep the satellites essentially 
"drag-free" by controlling radial, along-track and cross-track thrusters to negate the 
effects due to drag and solar radiation pressure, so that high sensitivity to the 
Earth's gravity field will be achieved [Keating etal, 1986]. 


It is desired that the groundtrack of each satellite repeats after a certain 
number of days that are not commensurate with the orbital period. The number of 
days required for the groundtrack repeat determines the resolution of global 
coverage. 

While in a frozen orbit the mean orbital ellipse will not change its eccentricity 
or orientation in space. This means that each satellite will fly over the same location 
on Earth at nearly the same altitude, and their respective perigee location will remain 
fixed over the north pole. 

The GRM launch date is currently proposed for sometime in the period 
1993-95. After launching one of NASA's Space Shuttles into polar orbit from the 
Western Space and Missile Center, the Remote Manipulator Unit (RMU) will 
deploy the two satellites from the Space Shuttle's cargo bay. At a nominal altitude 
of 275 km, both satellites will be deployed in the same orbital plane but separated 
by 50 km. Then after several Hohman transfers and daily checkouts, the satellites 
are lowered to the 160 km altitude with a separation of 150 km to begin mission 
operations 168 hours after deployment. 

The orbits of the spacecraft will be determined by tracking observations 
provided by the Defense Mapping Agency (DMA) with NASA's Tracking Data 
Relay Satellite System (TDRSS), the Global Positioning System, the TRANET 
ground based tracking network, or a subset of these techniques. 

At the start of the mission with 1400 kg of hydrazine fuel aboard to power 
the DISCOS and attitude control thrusters, the leading satellite would have a mass 
of 2734 kg while the trailing satellite has a mass of 2517 kg. This mass difference 
is due to a boom that carries the magnetometers, and star cameras which maintain 
the spacecraft’s orientation in space that are included on the lead spacecraft. Each 
satellite has a frontal area of 1.06 m 2 and an expected drag coefficient of 3.5. The 


external and internal configurations of the lead GRM satellite are shown 
respectively in Figures 1.2 and 1.3. 

1.2 The DISCOS System 

The atmosphere at 160 km is dense enough to cause rapid orbit decay. If 
drag effects are not counteracted, each satellite's lifetime would be less than 3 days 
[Holmes et al, 1987]. Since the density of Earth's atmosphere is directly related to 
the current solar activities such as the amount of solar flux received by the Earth, it 
can not be easily predicted. Furthermore, since drag and solar radiation pressure 
would corrupt the science data, a DISCOS system will be used aboard the low 
altitude spacecraft to ensure that these disturbances are kept to a minimum. This can 
be achieved by firing thrusters to contain an inner satellite "proof mass" within a . 
cavity at the satellite's centroid. 

Eventually, the outer portion of the satellite drops in altitude because of drag 
and moves ahead of the proof mass, thus causing the proof mass to move aft with 
respect to the center of mass of the outer satellite. Upon the proof mass's approach 
to a deadband threshold, the corresponding thrusters fire on the outer satellite until 
the proof mass has a predetermined velocity with respect to the satellite's mass 
center. Relative to the outer satellite, the proof mass then moves foward until drag 
causes the relative velocity of the proof mass to become zero and then reverses 
direction. The proof mass then approaches the deadband threshold again and a 
control limit cycle is maintained. This steady-state condition involves relatively 
short firing pulses and relatively long non-thrusting coast arcs. The details of this 
DISCOS system are presented later in Chapter 4. 



1.3 T he -SareU it £ztQTSate Hi t e Tr acking Sy st em 


Line-of-sight relative range-rates between the two satellite are provided by 
the Satellite- to-Satellite Tracking System. With this system, each satellite's 
transmitter continuously transmits millimeter radio signals at 42 GHz and 91 GHz to 
each other. These frequencies were chosen to keep errors due to ionospheric effects 
below the limit necessary to satisfy the accuracy requirements [Keating et al, 
1986]. 

As the pair of sateelites approach and pass a gravity anomaly, variations in 
the motion of the satellites are first induced on the lead satellite and then on the 
trailing satellite. As a result, the relative velocity between the satellites varies, from - 
1.0 m/s to + 1.0 m/s. The variation in velocity, referred to as the relative range-rate, 
is an indirect measure of the gravity anomalies. The range variations are measured 
by the Doppler-shift of the frequency of the received signal when compared to a 
reference signal provided by a 5 MHz oscillator [Keating et al, 1986]. 

1.4 Requirements of GRM 

The Satellite-to-Satellite Tracking system is required to measure the relative 
velocities between the antennae on each spacecraft to better than 1 (J.m / sec. This 
requirement is essential for the recovery of the gravity field, so that at an altitude of 
160 km the data set would provide resolutions of at least 1-2 mgal (1 mgal = 10 -5 
m/s 2 ) in gravity anomalies, and 5-10 cm in geoid height with a spatial resolution of 
100 km These resolutions are recommended by national scientific communities for 
the study of geological features such as sedimentary basins, batholiths, mountain 
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ranges, subduction zones and shields [Keating et al, 1986]. Also, the magnetic 
field will be resolved to 1 nT (nanotesla) at a 100 km spatial resolution [Keating et 
al, 1986]. The orbit determination accuracy (3a) requirements for the gravity 
mission are 100 meters in the radial direction and 300 meters in the along- and 
cross-track directions, whereas the magnetic mission requires an accuracy of 60 m 
in the radial direction, and 100 m in the along- and cross-track directions [Keating, 
etal, 1986], 

1 .5 Organization 

The research presented in this work covers details of the high precision and 
nominal orbit simulations which are discussed in Chapter 2. Chapter 2 also 
describes the mathematical formulae, software and hardware used for the generation 
of each simulation. The results of this study are shown in Chapter 3. The DISCOS 
instrument is discussed along with its history and implementation with GRM in 
Chapter 4. Also explained in Chapter 4 is the Drag-free Thruster Control Algorithm 
that contains all the control logic to operate DISCOS along with all the mathematical 
formulae, force models, software, and hardware that are used to simulate it. The 
results of this study are also presented in Chapter 4. The conclusions of this report 
are discussed in Chapter 5. 
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Figure 1.2 The External Configuration of the Lead GRM Satellite (Satellite 1) 

[Keating et al, 1986]. 



ORIGINAL PAGE IS 

of poor quality: 



— 

— 

— 





- — 

=a 





L 

zzi; 

__ 

— 

= 

— 

— 


— 

E 

E 

— 


— 

E 

— 

— 

— 

— 

— 

_ 






Figure 1.3 The Internal Configuration of the Lead GRM Satellite 

[Keating etal , 1986]. 





CHAPTER n 


SIMULATION OF GRM 

Two orbit simulations were performed for GRM. The first simulation 
consisted of a high precision representation of the actual orbits of the proof masses 
inside the cavities of the two low altitude satellites described in Chapter 1. The 
second simulation consisted of nominal orbits of the two 'inner' satellites 
determined by least squares approximation of the actual orbits. 

The high precision 'true' orbit was generated using an Ohio State University 
gravity field complete to degree and order 360 [Rapp, 1986], referred to as the 
"OSU86F" field. The ephemerides of the two low altitude satellites were generated 
at 4 second intervals for 32 solar days. Satellite-to-satellite line-of-sight range-rates 
measured by each spacecraft's integrated one-way Doppler antenna are also 
included in this simulation. The initial conditions of each satellite orbit were 
especially designed so that their groundtracks repeat after 32 sidereal days [White, 
1987]. This 32 sidereal day repeating groundtrack results in exactly 525 orbital 
revolutions for each satellite. Therefore, the orbital period is not commensurate with 
the repeat period so that high resolution coverage of Earth's gravity field can be 
achieved (1° by 1°). To concentrate on the effects due to high degree and order 
terms in Earth's gravitaional field, the gravitational effects due to the Sun or Moon 
were not included in the force model nor were effects due to kinematic or temporal 
forces. To keep computer costs low and maintain extreme accuracy, the KSG 
[Lundberg, 1984] fixed-mesh multi-step integrator was used along with the Encke 
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formulation of the differential equations. Furthermore, to avoid singularities at the 
Earth's poles, Pines' [1973] formulae that uses direction cosines rather than 
spherical coordinates were used to evaluate the geopotential. 

The nominal orbit includes the ephemerides of the two proof masses estimated 
by performing a least squares fit to the true ephemerides plus the estimated range- 
rates between the two 'inner' satellites. The nominal orbit was generated using a 
modified Goddard Earth Model 10B (GEM10B) gravity field that includes 
coefficients to degree and order 36 with the zonals J2 and J3 and the first two pairs 
of resonant coefficients at orders 16, 17, 33, 49, and 82 adjusted [White, 1987]. 
The initial conditions of the nominal orbits were especially designed for this 
simulation to maintain low radial, transverse and normal (RTN) differences 
between the true and nominal orbits (residuals), but its 32 sidereal day mission was 
not required to repeat [White, 1987]. The University of Texas Orbit Processor, 
UTOPIA, along with the Encke formulation of the differential equations and Pines 
[1973] formulation of the geopotential were used to generate the nominal orbit 
simulation. 

This chapter contains the mathematical formulae that model the motion of each 
GRM satellite in orbit about a non-spherical Earth. Section 2.1 presents the 
equations of motion that describe the orbit of each satellite and the force model that 
perturbs the satellites' motion. Pines [1973] formulation of the Geopotential is 
explained in Section 2.1.1. Increased numerical accuracy is obtained through the 
use of the Encke integration method which is presented in Section 2.1.2. The initial 
conditons for both the true and nominal simulations are discussed in Section 2.2. 
Section 2.3 lists the parameters that describe the Encke reference orbits. The gravity 
models are discussed in Section 2.4. The Integrated One-Way Doppler Range-Rate 



computer algorithm that simulates the Satellite-To-Satellite Tracking system is 
presented in Section 2.5. The equations of motion are integrated efficiently, saving 
computer costs with the use of the KSGFS numerical integration algorithm 
described in Section 2.6. Finally, the specifications of the Cray X-MP/24 
supercomputer used to numerically integrate the satellites’ orbits is described in 
Section 2.7. 


2.1 The Equations of Motion 


The orbit of each satellite is modeled by Newton’s Law of Gravitation which 
states that the force due to mass attraction between two point masses M and m 
separated by a distance r is given by the product of the masses, divided by the 
square of the distance and multiplied by a gravitational constant, G, or 


F = 


GMm 


( 2 . 1 ) 


By cancelling the satellite's mass (m) from both sides of Equation (2.1) and 
adding a perturbative force, the equations governing the motion of each satellite are 


r = — + f (2.2) 

r 

where 

r = the position vector of the satellite as referenced to a nonrotating, 
geocentric coordinate system, such as defined by a mean epoch of 
2000, 

|i = the gravitational parameter of Earth, 
f = the perturbing force due to nonsphericity of Earth. 



The proof mass will be susceptible to only conservative gravitational forces. 
Only the gravitational forces due to the Earth's static potential were examined. The 
forces due to third-body effects from the sun or moon, temporal effects from land 
and ocean tides, kinematic effects from precession, nutation, polar motion and UT1 
variations were not included in this model because they would complicate the 
gravity evaluation techniques. Nonconservative forces such as drag and solar 
radiation pressure, likewise, will not effect the proof mass because of the DISCOS 
system. The gravitational force is derived from the gradient of the Earth's potential, 
that is, 

f= VU (2.3) 

where the geopotential, U, is expressed in terms of spherical harmonics as: 

n=l 

0° OO a n 

+ -^XX(t) Pnm^^tC^cosmX + S^sinmX] (2.4) 

n=lm=l 

where r, <J>, X are the radius, latitude, and longitude of the satellite’s proof mass, 
modelled as a point mass with respect to the center of Earth, and p are the mean 
equatorial radius and the gravitational parameter of the Earth respectively, P^, are 
the associated Legendre functions, and S^, are the constant geopotential 
coefficients, and n and m are the degree and order of the geopotential field. 

2J .I Pines Formulation of the Geopotential 


Since the GRM mission specifies a polar orbit for global coverage of Earth, a 



singularity problem exists at the poles when the gradient of the spherical harmonic 
geopotential representation in Equation 3.4 is calculated. By using direction cosines 
rather than spherical coordinates. Pines [1973] reformulated Equation 2.4 so that 
the singularities no longer exist. The geopotential can now be used for polar orbits 
with the form being 

vva s n 

U + A ( u ) C 

j* j* f ' n N ' n,o 
n=l 

oo oo 

yr ' a n 

+ a (u)[ C r (s,t) + S i (s,t)] (2.5) 

r *■'« « ■■■'< v r ' nm v / 1 n,m m' ’ ' n,m m v ’ /J v ' 

n=lm=l 

where s, t, and u are the direction cosines, 



A^ m are the "derived" Legendre functions, and r m (s,t) and i m (s,t) are the real and 
imaginary parts of (s + j t) m . Recursion algorithms used for the calculation of the 
"derived" Legendre functions in Pines' formulation on high speed computers is 
discussed by Lundberg and Schutz [1987]. 

2.1.2 The Encke Formulation 

A technique to perserve numerical accuracy during integration is provided by 
the Encke integration method. This technique uses an analytic reference orbit which 
closely matches the true orbit. The Encke vector then is described as the difference 
between the satellite's true trajectory and the reference orbit. By numerically 
integrating the Encke vector rather than the satellite's true state several significant 
digits can be perserved [Lundberg, 19S5]. 
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The Encke vector denoted by £ is equal to the difference between the reference 
orbit, r s and the true orbit position, r as shown in Figure (2.1). The Encke vector 
with its first and second time derivatives are then 



( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 


The reference orbit for this work is modelled as a secularly processing ellipse 
and can be written in differential form as 




mj - 

r = - — i + f (2.9) 

s 3 s v / 

r 

s 

where r s is the position vector to the reference orbit in a non-rotating, geocentric 
coordinate system as in Equation (2.2) and f s is the perturbing accelerations that 
result in a secularly precessing ellipse. The substitution of Equations (2.2) and 
(2.9) into Equation (2.8), results in the Encke formulation. 


llr llr _ 
e = - + —1 +f-f 

3 3 * 

r r s 


( 2 . 10 ) 


Extreme care was taken to define the reference orbit for these simulations, 
since possible problems may result in this formulation if the reference orbit nearly 
equals the true orbit or if it does not remain close to the true orbit. With the first 
problem, a signifcant loss of numerical accuracy may result from the subtraction of 


nearly equal terms. However, this problem can be overcome by arranging the 
computations for the first two terms of Equation 2.10 in an appreciate form. To 
achieve increased accuracy, the Encke vector, £, must remain significantly smaller 



than either the true or reference position vectors which means that the reference 
orbit must not deviate substantially from the true orbit. The choice for the reference 
orbit used in these simulations was provided by Lundberg [1986]. 

2.2 Initial Conditions 

The initial conditions (x = 0.) for the true and nominal ephemerides were 
provided by White [1987]. The true orbit's initial conditions were especially 
designed so that the final state after exactly 32 sidereal days will be very nearly the 
same as the initial conditions. They were determined using the University of Texas 
Orbit Processor (UTOPIA) and are given below in Table 2.1. The initial conditions 
for the nominal orbits shown in Table 2.2 were also estimated using UTOPIA, by 
fitting the emphemerides of the true orbits in a least squares sense with a different 
model to simulate an actual situation. 

Table 2.1 Initial Conditions for the 'True 1 Orbit 

Satellite 1: 

X} = 26116184162 m, Y t = -150104.5682242 m, = 6515224.696995 m, 

V X j = -4.81851974 x 10* 2 m/s, V Yl = -7816.577574349 m/s, V Zl = -179.5770526472 m/s, 

Satellite 2: 

X 2 = 262.89992177 m, Y 2 = 149884.9023112 m, Z 2 = 6515227.869697 m, 

V X 2 = -4.76637746 x 10‘ 2 m/s, Vy 2 = -7816.587219218 m/s, Vz 2 = 179.4189749253 m/s, 


Table 2.2 Initial Conditions for the Nominal Orbits 


Satellite 1: 

Xj = 251.82011984 m, Y x = -150205.8333573 m, Z t = 6515206.75655362 m, 

V X j = -8.617240919 x 10‘ 2 m/s, V Yl = -7816.5813004031 m/s, V Zj = -179.65789276711 m/s, 

Satellite 2: 

X 2 = 254.34316160 m, Y 2 = 149823.0882622 m, Z 2 = 6515206.44534594 m, 

Vx 2 = -8.650294014 x 10' 2 m/s, Vy 2 = -7816.6011900852 m Is, \z 2 = 179.39406524579 m/s 

2.3 Encke Reference Orbits 

The numerical integration of both the true and nominal orbit simulations used 
the same Encke reference orbits that are described in Section 2.1.2. The orbital 
parameters used to generate these reference orbits are listed in Table 2.3. 


Table 2.3 Encke Reference Orbit Mean Orbital Elements 



Satellite 1 

Satellite 2 

Semimajor Axis, a (m) 

6523600.811305 

6523599.627289 

Inclination, i (deg) 

90 

90 

Longitude of Ascending Node, fl (deg) 

90 

90 

Argument of Periapse, m (deg) 

0.0 

0.0 

Eccentricity, e 

0.0 

0.0 

Mean Anomaly, M (rad) 

1.593311064614 

1.547312542033 

Node Rate, £1 (rad/sec) 

0.0 

0.0 

Argument of Perigee Rate, to (rad/sec) 

0.0 

0.0 

* 

Modified Mean Motion, M (rad/sec) 

1.1963632130262 x 10* 3 

1.1963633652585 x 10’ 3 


2.4 The Gravity Model 


The gravity model, OSU86F, used in the GRM high precision simulation was 
created at Ohio State University [Rapp & Cruz, 1986]. It contains a set of constant 
geopotential coefficients complete to degree and order 360. The gravitational 
parameters used in this model are given in Table 2.4. 

Table 2.4 Gravity Model for the True Orbits 


Initial Greenwich Hour Angle, 

a = 

100.3399460 deg, 

Earth's Rotational Velocity, 

“e = 

7.2921158553066 x 10' s rad/s 


= 

[15.04106864 deg/hr], 

Earth's Gravitational Parameter, 

p = 

3.986004404 x 10 14 m 3 /s 2 , 

Earth's Mean Radius, 

a e = 

6378137.0 m, 

Maximum Degree of Geopotential: 


360, 

Maximum Order of Geopotential: 


360 


The Goddard Earth Model 10B (GEM10B) used for the nominal orbit was 
created at the Goddard Space Center and is discussed by Lerch et al [1981]. 
Coefficients up to degree and order 36 plus some additional coefficients (up to 
order 82) that produce resonance effects on the GRM orbit were included also. The 
gravitational parameters for this model are: 


Table 2.5 Gravitational Model for the Nominal Orbits 


Greenwich Hour Angle, 

Earth's Rotational Velocity, 

Earth's Gravitational Parameter, 
Earth's Mean Radius, 
Geopotential: 


a = 100.3399460 deg, 

( 0 e H 7.2921158553066 x 10' 5 rad/s 

= [15.04106864 deg/hr], 

|i = 3.9860064 x 10 14 m 3 /s 2 , 

a e = 6378145.0 m, 

G10B complete to degree and order 36 
with adjusted resonant coefficients of 
orders 16, 17, 33, 49 and 82 


2,5 The One-Way Doppler Range-Rate Algorithm 


Gravity anomalies will be indirectly measured by the Satellite-to-Satellite 
Tracking System (SST) that uses an integrated one-way Doppler system to calculate 
range-rates between the proof masses of the two satellites. Simultaneously, while 
the lead satellite's (A) narrow-beam antenna transmits two signals at 42 and 91 
GHz directed towards the trailing satellite (B), the trailing satellite’s antenna 
transmits two signals at the same frequencies towards the leading satellite. These 
two frequencies, that are broadcast continuously, were chosen to filter out 
ionospheric effects. When the line-of-sight distance between the satellites varies, 
the frequency of the signals are doppler shifted. These doppler-shifts are measured 
by comparing the signals to a 5 GHz reference signal. Satellite-to-satellite range- 
rates are then computed from the doppler-shifts. Finally, gravity anomalies can be 
evaluated using these range-rates with various techniques. 
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The one-way doppler algorithm simulates the doppler measurement by first 
calculating the time of flight (TOF=Ar) of each signal of the same frequency. If 
both satellites receive each other's signal at q then the trailing satellite (B) 
transmitted a signal to the lead satellite at time q - At Bi and the leading (A) 
transmitted a signal to the trailing satellite at time q - At Ai . The path length or the 
distance between the satellites is calculated by subtracting the position vector of the 
transmitting satellite at time q - At i from the position vector of the receiving 
satellite at time q. This requires interpolating the states of each satellite at the time 
of transmission. For each satellite, the range or distance , p } , that their signals have 
travelled must be iterated until the following equations hold true, 

PBt = I f B<' I = c < At Ai> < Z12 > 

where At is the TOF of the corresponding signal and c is the speed of light in a 
vacuum. An average of each satellite's simultaneously measured distance is then 
taken so that the path length of the signal that arrives at time q is 

f*Ai + PBi 

Pi = j (2-13) 

Finally, the integrated range-rate measurement at qis determined by 
subtracting the previously measured range from the current range value and 
dividing by the interval between measurements, h, which for this simulation was 
taken to be four seconds, i.e.. 



( 2 . 14 ) 


The analytical equation for the instantaneous range-rate is 
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PlNST “ 


PlNST ' PlN ST 
PlNST 


(2.15) 


where the instantaneous range vector is the difference between the position vector 
of the leading (r A ) and trailing satellites (r B ), 


Pinst r A ‘ r B (2.16) 

PlNST * s r ^ e magnitude of the instantaneous range vector, and is the time 
derivative of Pinst- 

Instantaneous range-rates from Equation (2.15) were calculated along with the 
integrated one-way doppler range-rates from Equation (2.14) and compared. These 
results are presented in Chapter 3. 


2.<L The KSGFS Numerical Integrator 

A time saving, but highly accurate numerical integration algorithm, KSGFS 
[Lundberg, 1985], was used for the creation of the ephemerides of all orbit 
simulations . The KSGFS package is a multi-step, fixed mesh, class II numerical 
integrator that uses a constant stepsize (h) with the general formulation (Adams' 
type) of the Predict-Evaluate-Correct-Evaluate (PECE) Algorithm. The accuracy of 
the KSGFS integrator is a function of the order of the integrator (NORD), the 
constant step size (h), and the starting convergence criterion (ALIM). The order of 
the integration, which represents the number of function evaluations that are stored 
within the back-difference table, was set at 10, the stepsize was set at 4 seconds, 
and the starting convergence criterion was set at 10* 16 . 


2.7 The Computational Facilities 


All computations within this study were carried out on the University of 
Texas System Center for High Performance Computing (CHPC) Cray X-MP/24 
supercomputer located at the Balcones Research Center in Austin. The CFT 1.15 
fortran compiler under the Cray Operating System (COS) was used to process all 
software. Using a 64 bit word with a 48 bit mantissa, the Cray X-MP/24 provides 
approximately 14 decimal digits of numerical accuracy in single precision arithmetic 
operations. 

Since the evaluation of the Earth's geopotential requires 90 percent of the 
CPU computational time, it is essential that its recursion algorithms are completely 
vectorizable. Studies of these algorithms were performed prior to the actual 
simulation study by Schutz et al [1987] at Cray's Mendota Heights computer 
facility on the CRAY X-MP/48 and on the Cray X-MP/24 at U.T.. From these 
tests, the 32 day high precision orbit simulation with a geopotential of degree and 
order 360 was estimated to require approximately 20 hours of Cray's CPU time for 
a numerical integration step size of 4 seconds. Actual simulation time costs are 
reviewed in Chapter 3. 
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CHAPTER m 


SIMULATIONS OF THE TRUE AND NOMINAL ORBITS 

This chapter presents the results of creating the true and nominal orbit 
simulations performed at the University of Texas at Austin. Section 3.1 describes 
the predicted computer costs (CPU time) involved in creating a high precision orbit 
simulation that will represent the true or actual GRM orbits. The actual computer 
cost is discussed in Section 3.2 along with the results of the true orbit simulation. 
The comparison of the final conditions with the initial conditions is also presented 
in Section 3.2 to illustrate how closely they fullfill the groudtrack repeat 
requirement 

3.1 Pre-evaluation Computer Cost Analysis 

Before submitting the costly High Precision Orbit Simulation, a computer cost 
analysis was performed on the software that was used generate it. It was found that 
evaluating the gravity accelerations by using the geopotential routine (GRAV) takes 
90 percent of the execution costs (CPU time) on the Cray X-MP/24. The cost 
analysis involved timing a number of executions with the GRAV routine using 
different geopotential files in both single and double precision. Three geopotential 
files were investigated in this study: OSU322 {Rapp, 1981], OSU300 and 
OSU86F {Rapp et al, 1986] from the Ohio State University. The goal was to find 
highly accurate but efficient software to evalutate the geopotential. The results are 
presented in Table 3.1. 
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Since the Cray Fortran compiler does not vectorize double precision 
arithmetic, the timings per function evaluation are rather large. The OSU86F 
geopotential field was chosen for the high precision orbit simulation because it was 
the most recent file published that contains higher degree and order coefficients 
based on actual observations of surface gravity. With a 4 second integration 
stepsize and 2 function evaluations per integration step, the 32 day simulation was 
expected to take nearly 20 hours of CPU time on the Cray X-MP/24. 

3.2 The True Orbit Simulation 

With a four second integration step size, the 32-day repeating groundtrack 
mission required 19.2 hours of CPU time for the computation of 1.38 million 
function evaluations on the Cray X-MP/24 supercomputer. The final state (position 
and velocity) of the two "inner" satellites (proof masses) after 32 sidereal days is 
presented in Table 3.2. The ephemerides of both satellites were recorded on 
magnetic tapes at 4 second intervals for the 32 solar days. Also recorded on 
magnetic tapes were the Encke vector along with its first and second derivatives 
which represent the position, velocity, and acceleration with respect to the reference 
orbit; the true position and velocity vectors; the acceleration vector due to the 
nonspherical geopotential perturbations; and the Satellite-to-Satellite Tracking 
range-rate measurements. The 32 solar day ephemeris includes 691,200 data 
points. 
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By comparing the initial longitude and latitude with the final longitude and 
latitude of each satellite, it was determined that the orbits repeat at the end of 32 
sidereal days to within 1.78 km for the lead satellite and 2.52 km for the trailing 
satellite. This means the lead satellite misses its repeat location by 0.233 seconds, 
and the trailing satellite misses its repeat location by 0.331 seconds. As shown in 
Table 3.3 which compares the initial and final conditions of each satellite, the 32 
sidereal day orbits resulted in very close repeats of their initial conditions, thus 
satisfying mission requirements. 


Table 3.3 Comparis on of the Initial and Final Conditions 



Satellite 1 

Satellite 1 

Initial Latitude: 

88.68881 deg 

88.69073 deg 

Final Latitude: 

88.67285 deg 

88.71340 deg 

Latitude Error: 

0.015960 deg 

0.022400 deg 

Longitude Error: 

0.000600 deg 

0.001640 deg 

Groundtrack Error: 

1.780 km 

2.5240 km 

Time Error: 

0.233 sec 

0.33072 sec 


3.3 The Nominal Or bit Simulation 

The nominal simulation was determined by performing a least squares fit to 
the true ephemerides described above using the University of Texas Orbit Processor 
(UTOPIA). The fit required the adjustments of the zonals J 2 and J 3 and the first two 
pairs of the resonant coefficients at orders 16, 17, 33, 49, and 82 in the GEM10B 
gravity field. The initial conditions of the nominal orbit simulation were estimated 


along with the above resonant coefficients by White [1987] and were presented in 
Chapter 2. The final conditions of the nominal orbit simulation at the end of 32 
sidereal days are presented in Table 3.4. 

It was desired to use a small geopotential field differing from the thruth model 
to simulate an actual mission. Such a field results in differences that exceed the 
mission requirements. Therefore, the nominal simulation was first estimated for the 
leading satellite using a 36 by 36 GEM 1 OB gravity field. This baseline trajectory 
was then differenced from the true trajectory of this satellite in the radial, along- 
track, and cross-track directions. These differences were rather large and they did 
not meet the mission’s specifications of 100 m in the radial and along-track 
directions and 300 m in the cross-track direction. It was determined that the 
spherical harmonic coefficients at orders 16, 17, 33, 49, and 82 were in resonance 
with the 160 km orbit since their periods are commensurate with an integer number 
of the satellite’s daily revolutions. These resonant coefficients were seen to cause 
periodic variations in the along-track direction. For example, order 82's resonant 
period was found to be 32 days with an amplitude over 800 meters [White, 1987]. 
Therefore, the first two pairs of the harmonic coefficients, and S^,, were 
estimated at these orders to reduce the residuals within the mission's specifications. 
However, these specifications were still not completely met so it was decided that 
the zonals J 2 and J 3 should also be estimated. Once the trajectory of the leading 
satellite was estimated with residuals to meet the mission specifications, the 
adjusted GEM10B gravity field was used to estimate satellite two’s trajectory. The 
details of this estimation process for the nominal orbit are explained by White 
[1987]. 






3.4 RTN Differences Between True and Nominal Orbits 


The true and nominal orbit simulations are compared by examining the radial 
(R), transverse (along-track, T), and normal (cross-track, N) components of the 
position vector difference along the entire trajectory of the 32 day orbit. These 
differences indicate the level of effort necessary to meet the orbit determination 
requirements of the mission. The prime concern of this study is the gravity mission 
which requires 3a accuracy of 100 meters in the radial and 300 meters for along- 
and cross-track components. The differences between the true and the nominal 
orbits must remain below these requirements. The radial component of the 
differences versus time for the lead and trailing satellite (satellite 1 and satellite 2) 
are shown in Figures 3.1-a and 3.1-b, respectively. The transverse difference 
versus time are shown in Figures 3.2-a and 3.2-b for satellite one and satellite two, 
respectively. The normal differences are shown in Figures 3.3-a and 3.3-b. These 
graphs were generated by sampling data at a rate of five points per orbit. 

As seen from these figures, the radial differences of both satellites remain 
between ±55.0 meters; the transverse residuals of both satellites remain between 
+200.0 and -250.0 meters; and the normal differences of both satellites remain 
between ±55.0 meters. Therefore, the nominal orbit simulation meets the 
requirements for the gravity mission, whereas the magnetic mission requirements 
(3 a) of 100 meters in the along-track direction are not satisfied. However this 
nominal orbit simulation will be used to improve gravity models so that eventually 
the magnetic mission's requirements can be met, or additional resonance terms can 
be estimated to further reduce the transverse differences. 



3.5 Range-Rate Measurements 


If the range-rate between the satellites could be instantaneously measured, 
then a plot of range-rate versus time for the 32 day true orbit simulation would 
have the characteristics shown in Figure 3.4. Range-rates during the beginning of 
the mission would range between ±0.8 m/s, then increase gradually (±0.9 m/s) near 
the end of 32 days. 

Actual range-rate measurements made by the SST's integrated one-way 
Doppler system along the true orbit would resemble those shown in Figures 3.5-a 
and 3.5-b for the reception of the signal by satellite 1 and satellite 2, respectively. 
These plots also exhibit the trend described in the instantaneous plot. These 
simulated measurements were made by the one-way Doppler range-rate algorithm 
discussed in Chapter 2. Figures 3.6-a and 3.6-b show the time of flight (minus 1.0 
millisecond) of the signal received respectively by satellite 1 and satellite 2 against 
time. Since the variations of the time of flight are on the order of 1.0 ps, one 
millisecond is subtracted from the actual time of flight in order to show these small 
differences in the time of flight. With the speed of light (c) approximately 3 x 10 8 
m/s, the corresponding distances between the satellites could also be represented by 
these figures. The average range-rates of these two measurements are presented in 
Figure 3.7, and the differences of the two measurements are presented in Figure 
3.8. The two satellites' antennae are not measuring the same range-rates since each 
satellite is accelerating and decelerating at different times and signals are not 
instantly received. Figure 3.8 shows the differences between the two satellite 
measurements to range between ±100 microns per second. 

All measurements calculated by the one-way Doppler algorithm were 
compared against the ideal instantaneous measurements. Figures 3.9-a and 3.9-b 


show the differences between the measurements calculated from the signal received 
using the one-way Doppler algorithm and the instantaneous measurement for 
satellite 1 and satellite 2, respectively. Figure 3.10 shows the difference between 
the instantaneous and the average of the two satellite range-rate 
measurements.These differences remained between +2.5 and -3.5 millimeters per 
second. 

The nominal orbit simulation range-rate measurements were then compared to 
those from the true orbit It is anticipated that the difference between the nominal 
and true measurements would be used to recover the true geopotential field. Figures 
3.1 1-a and 3.1 1-b show the differences between the nominal and the true range-rate 
measurements against time as measured respectively by the signal received by 
satellite 1 and satellite 2. These differences remain below ±35 millimeters per 
second. 

Each of the above figures were generated by sampling the data at a rate of five 
points per orbit. The sinusoidal patterns seen on these graphs and the above RTN 
graphs are usually the results of this sampling rate and are generally not the true 
indication of the actual data. 

It was found that the iterating procedure of the one-way Doppler algorithm to 
determine the time-of-flight of the transmitted signal of satellite two was unable to 
converge occasionally (15 points out of a total of 691,200 measurements). The 
computer code and algorithm was extensively evaluated to clear this problem, but 
the cause was not determined. Possibly, it could be due to round-off error when 
interpolating for the states of each satellite. The measurements calculated during 
these non-convergences did not meet the 1 |im/s accuracy required, and were 
ignored. 
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Figure 3.5-a Integrated One-Way Doppler measurements for signal received by 
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3.5-b Integrated One-Way Doppler measurements for signal received by 
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Figure 3.6-a Time of flight of signal received by satellite one. 
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Figure 3.6-b Time of Flight of signal received by satellite two. 
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Figure 3.7 Average of both satellite's Integrated One-Way Doppler measurements 
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Figure 3.8 The difference between the integrated one-way doppler range-rate 
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Figure 3.9-a The difference between satellite one and instantaneous range-rates. 
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Figure 3.9-b The difference between satellite two and instantaneous range-rates. 


0.003 


1 

u 

04 



1 I 1 1 I ■’ 1 I ■ ■ I 1 ' I » ' » I ' 1 I 1 1 I ' ' I 1 ' 

6 9 12 15 18 21 24 27 30 33 


TIME (DAYS) 


DIFFERENCE (M/S) 


! INAL PAGE IS 

ur POOS QUALITY 


Figure 3. 10 The difference between instantaneous and average range-rates. 
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Figure 3.1 1-a Differences between true & nominal Integrated One-Way 
Doppler measurements for satellite one. 



0.040 


0.030 


0.020 

5? 


1 

0.010 

w 


o 

z 

0.000 

§ 


1 

- 0.010 

Q 

- 0.020 


- 0.030 


- 0.040 



I 1 ■ i ■ i ■ ■ i ■ ■ i ■ ■ i ■ ■ i ■ ■ i ■ ■ i ■ ■ i ■ i i i i | 

0 3 6 9 12 15 18 21 24 27 30 33 

TIME (DAYS) 


Figure 3.11-b Differences between true & nominal Integrated One-Way 
Doppler measurements for satellite two. 
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CHAPTER IV 


SIMULATION OF THE DRAG COMPENSATION SYSTEM 

This Chapter describes the mathematical formulae and control logic for the 
design of a Disturbance Compensation System (DISCOS) thruster control 
algorithm. The DISCOS is designed to offset effects due to drag and radiation 
pressure on the satellite, but only drag will be considered here. 

A description of the DISCOS mechanism along with its history is discussed in 
Section 4.1. Section 4.2 describes the control logic of the computational algorithm 
used for this simulation. The equations of motion for this analysis are derived in 
Section 4.3 and Section 4.4 descibes the force model used in this study. The 
numerical integration of the equations of motion of the proof mass and outer 
satellite is discussed in Section 4.5. The control laws, system logic and parameters 
for on-off thrusting switches is detailed in Section 4.6. Section 4.7 describes the 
ball-centering attitude control law that is coupled with the DISCOS translational 
control system so that small lift and sideslip forces can help reduce fuel 
expenditure. Finally, the results of this investigation are presented and the amount 
of fuel required to fire the thrusters for drag compensation during the six month 
mission is calculated in Section 4.8. Also, fuel expenditures are compared with 
previously published reports of Keating et al [1986] and Ray & Jenkins [1981] in 
Section 4.8. 
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4A The Disturbance Compensation System 


As early as 1962, the feasibility of a mechanism that freed a satellite from 
atmospheric effects was first investigated by B. Lange, [1962]. Later, a "drag- 
free" mechanism, the Disturbance Compensation System (DISCOS), was 
developed in a joint effort by The Johns Hopkins University Applied Physics 
Laboratory and the Standford University Guidance and Control Laboratory 
[1974] for the 1972 U.S. Navy's TRIAD satellite in an effort to improve navigation 
by satellite. 

In the DISCOS design, a small satellite (the proof mass) is completely 
enclosed inside a cavity at the centroid of a larger outer satellite. This inner proof 
mass is then shielded from all nonconservative forces such as drag and radiation 
pressure. As the outer satellite is perturbed by these forces, its orbit deviates from 
the proof mass's orbit. When its deviation from the proof mass reaches a 
predetermined limit (the deadband region), appropiate thrusters fire from the outer 
satellite's propulsion system to keep its orbit coincident with the inner satellite's 
orbit and effectively negates the effects of nonconservative forces. This design was 
observed to maintain the proof mass of the TRIAD satellite within specified 
boundaries at an altitude of 700 km for 18 months using only three pounds of cold 
propellant [JHU APL et al, 1974]. 

While in orbit, the spherical proof mass orbits without contact with the outer 
satellite's spherical cavity, and its position relative to the outer satellite is monitored 
by a set of capacitance bridges (2 plates oppositely aligned on each of the 3 
coordinate axes). External forces displace the outer satellite's electrical center of the 
capacitive cavity from the mass center of the proof mass, so that its capacity 


changes as a function of its position relative to the proof mass's center of mass. 

The TRIAD proved that a satellite equipped with the DISCOS mechanism can 
be essentially drag-free. Because a comprehensive gravity study such as GRM 
would require a low altitude orbit for extreme sensitivity to the geopotential, 
unpredictable, orbit decaying drag effects need to be constantly removed by an 
onboard control system. A review of the TRIAD DISCOS design indicated that its 
technology is adaptable to the GRM mission. Studies have shown that 0. 1 p.m/s 
minimum sensivity to the velocity of the proof mass with respect to the satellite’s 
centroid can be achieved by the modification of the TRIAD DISCOS sensor 
[Keating et al , 1986]. 

The proposed GRM DISCOS mechanism is comprised of a 160 mm diameter 
spherical beryllium oxide housing (cavity) that includes a capacitive three-axis 
orthogonal coordinate position sensor and associated electronics while enclosing a 4 
kg, 140 mm spherical aluminum proof mass and is located at the mass center of 
each satellite. Each entire DISCOS instrument, shown in Figure 4.1, is estimated to 
weigh 15 kg [Keating etal, 1986]. 

The outputs of the DISCOS sensor are used by the Guidance and Control 
System which controls the thrusting on all three body-fixed orthogonal axes. Each 
spacecraft uses eight along-track thrusters (4 foreward and 4 aft) that are rated at 4 
newtons apiece; eight 1 N thrusters (4 port and 4 starboard) for yaw momentum 
unloading and cross-track drag compensation; eight 1 N thrusters (4 top and 4 
bottom) for pitch momentum unloading and vertical translation control; and eight 1 
N thrusters (4 oppositely aligned) for roll momentum unloading. All thrusters are 
modelled as "idealized", instantaneous on-off force generators; and are fired in 


pairs. Approximately 1400 kg of liquid hydrazine provided by fuel tanks placed 
equidistant from the proof mass to minimize mass attraction will be used to fuel 
each spacecraft’s propulsion system. This gives a total impuse of 2.8 million 
newton-seconds of thrust to counteract drag and control each satellite's attitude 
[Ray, & Jenkins , 1981]. 

Since the Satellite-to-Satellite Tracking antennae are attached to the outer 
satellites, their precise displacement must be determined so that the relative velocity 
between the phase centers of the antennae can be determined to correct the SST 
data. As shown in Figure 4.2, the true range-rate (sj) is then 

&T = s + 5 x i + 6 x2 (4.1) 

• • 

where s is the measured range-rate provided by the SST and 5 xl and S x2 , 
respectively represent relative velocities of the leading and trailing spacecraft's SST 
antenna with respect to the proof mass [Keating et al, 1986]. 

According to Keating et al [1986], "the performance requirements of the 
DISCOS for the GRM mission are: 

1. ) Proof mass-to-spacecraft relative velocity knowledge within 0.1 jim/s 
(la) with a 4 second averaging time. 

2. ) Three-axis disturbance compensation within 10' 9 g, rms, over the 
frequency range of 0.001 to 0.5 rad/s." 

In order to meet the 10 -9 g compensation requirement, the proof mass (ball) 
must fly freely within the satellite's cavity without hitting its walls, or in other 
words, the outer satellite must follow the proof mass without touching it. 
Satellite/ball interactions due to mass attraction, electric and magnetic fields, 
residual gas pressure (inside the cavity) and radiation pressure are some of the 



possible sources of error. The 1972 flight of the U.S. Navy TRIAD satellite proved 
that these disturbances can be kept below 10' 11 g. It has been determined that 
motions caused by thruster firings will directly corrupt science data (Equation 4.1) 
unless it is detemined to a 0.1 |i.m/s accuracy. Ball center of mass offset, ball shape, 
pick-off noise, thermal deflections and structural flexure are some of the errors 
associated with the proof mass velocity measurements provided by the capacitive 
sensor. Also, since the capacitors (one pair on each axis) exhibit nonlinearities 
outside a small linear range (~±1.0 mm), the deadband limit is fixed at these 
boundaries. Therefore, the maximum distance the outershell can deviate from the 
proof mass is 1 mm (furthest ball excursion in all directions) before the 
corresponding thrusters turn on [Keating et al, 1986]. 

4.2 The Design of the Drag-Free Thruster Control Algorithm 

The primary objective of the thrusting control algorithm is to compensate for 
the drag effects on the outer satellite while minimizing fuel consumption. It involves 
controlling the outer satellite's thrusters so that it precisely follows the drag-free 
orbit of the proof mass within required limits. 

Both the along-track and cross-track thruster control models involve turning 
the appropiate thrusters on and off at the precise moments to keep the outer 
satellite's center of mass from deviating no more than 1 mm from that of the proof 
mass. The 160 mm diameter cavity allows the 140 mm diameter proof mass to drift 
within a 10 mm margin before it touches the wall. In general, the capacitance is a 
nonlinear function of proof mass displacement in the cavity. However, since the 
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compacitive sensors behave linearly within a ±1 mm region, the deadband zone (the 
thruster-on switch) is set at these limits. These on-off switches cause a contactor 
control limit cycle to occur in either the along-track (Figure 4.3) or cross-track 
(Figure 4.4) phase-plane diagrams. The thruster-on portions of the along-track limit 
cycle is represented by a parabolic arcs at the deadband regions in the phase-plane 
diagram as shown in Figure 4.3. The thruster-off portion is represented by a "coast 
arc" parabola in the phase-plane. The thruster-on/off switches occur at the 
intersection of these two parabolas. There is no steady-state limit cycle for the 
cross-track control law since there is no consistently predictable force in this 
direction. This law imparts a small inward velocity proportional to the proof mass 
distance from the deadband zone. 

As drag affects the outer portion of the satellite, the proof mass moves 
forward with respect to the center of mass of the outer satellite with increasing 
speed. Eventually the proofmass approaches the forward deadband region inside 
the cavity which is represented by the parabolic thrust arc in Figure 4.3. The aft 
thrusters turn on for approximately 60 milliseconds causing the outer satellite to 
move ahead of the proof mass. Then drag slows the velocity of the proof mass with 
repect to the satellite's centroid to zero and the ball reverses direction and speeds 
back to the forward deadband region as shown by the parabolic coast arc in Figure 
4.4. 

The frequency and amplitude of the along-track limit cycle was chosen to 
enhance the accuracy of the post-flight ball velocity reconstruction. Due to 
uncertainties, satellite-to-satellite range-rate measurements taken during thruster-on 
times are not to be used in the gravity recovery. Therefore, a repeatable limit cycle 
is needed for filtering out the data recorded during these times. However, the 
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variation of drag along the satellite's orbit cause difficulties in achieving a precisely 
repeatable limit cycle. Also it is desirable to maintain a small amplitude limit cycle to 
remain within the linear range of the sensor, but with a small amplitude cycle, the 
thrusters fire more frequently causing a increase in fuel consumption and also 
affecting thruster reliability. Therefore, "the desired steady-state condition is a 
sequence of relatively short firing pulses and relatively long coast arcs," [Ray and 
Jenkins , 1982]. 

Unwanted disturbance forces such as electrical charge and self-gravity on the 
proof mass are a important sources of error and the control law should minimize 
their effect on science data. Therefore, the average value of the proof mass 
excursions must be kept as close to zero as possible. According to Ray and 
Jenkins, [1982], "the average position on the parabolic arcs will be zero for x max - 
2 x^, which has been chosen for the nominal along-track limit cycle on GRM," 
where Xnj^ is the maximum along-track excursion of the proof mass with respect to 
the outer mass (apex of thrust arc parabola) and x mill is the minimum along-track 
excursion (apex of the coast arc parabola). So with these considerations, "the GRM 
limit cycle has been tentatively set at x^ = +1.0 mm and = -0.5 mm." 

Since the spacecraft's attitude affects the translational dynamics of the proof 
mass relative to the outer spacecraft through both aerodynamic forces and kinematic 
pseudo forces (centrifugal and corriolis) and the attitude thrusters will produce a 
torque and a force when fired, the attitude and translational control laws must be 
coupled [Ray and Jenkins, 1982]. The ball-centering attitude control system 
explained in Section 4.7 is an example of this coupling. 


52 


4.3 The Equations of Relative Motion 


Figure 4.5 shows the relative position of the inner satellite (proof mass) with 
respect to the outer satellite's center of mass. Both inner and outer satellites are 
modelled as point masses, however it is assumed that the mass attraction between 
the two is negligible. In a mean of 2000, equatorial, geocentric, inertial coordinate 
frame the position of the proof mass with respect to the outer satellite (outer mass) 
is 


e 


= r 

pm 


r 

om 


(4.2) 


where is the position vector of the proof mass and r om is the position vector of 
the outer mass as referenced to the inertial coordinate system. 

Taking the time derivative twice gives the proof mass's acceleration with 
respect to the outer mass as 



where the dynamical motion induced on the proof mass state is purely gravitational 
and is represented by the two-body acceleration plus accelerations due to a 
nonspherical Earth geopotential of degree and order 5, 
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The dynamical motion of the outer mass includes the two-body and 
nonspherical accelerations, but is perturbed by drag forces that are eventually 
counteracted by the thrusters. Therefore, the equations of motion of the outer mass 
is 
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r 

om 


= -|i- 


om 


+ f + f 


drag 


f + f 
1 thrust a BC 


(4.5) 


In the above equations, |l is the Earth’s gravitaional parameter, f gpm and f gom 
respectively represent the gradient of the geopotential induced on the proof mass 
and the outer mass due to Earth's nonsphericity. Drag and thrust are respectively 
represented as f drag and f^ust- The term f BC is a small lift force in the radial 
direction or a sideslip force in the cross-track direction produced by pitching or 
yawing the outer mass at small angles with respect to the velocity vector, so that 
they will help center the proof mass in these directions. This ball-centering attitude 
control law is explained further in Section 4.7. 

Finally, by substituting Equation 4.3 and Equation 4.4 into Equation 4.2, the 
equations of motion that govern the proof mass's state with repect to the outer mass 
are produced. 



(4.6) 


The state vectors of the proof mass with respect to the inertial frame and with 
respect to the centroid of the satellite are propagated by numerically integrating 
Equations 4.4 and 4.6 as opposed to Equations 4.4 and 4.5. The Encke integration 
can be used to provide greater accuracy in the solution for the relative motion than 
the differencing of the numerical solution of the absolute positions. 
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4.4 The Force Model 

To keep this model as simple as possible, forces due to three-body, temporal, 
and kinematic effects were not included in this analysis. Disturbance forces such as 
electrical charge on the proof mass and ball-satellite gravitational attraction are also 
not included in this analysis. 

4.4.1 Geopotential Forces 

The force induced on either the proof mass or the outer mass from the Earth's 
nonsphericity is respectively derived by the gradient of the geopotential, 

f , = TO „ < 4 - 7 > 

pm 

f, = ^on, «.8) 

om 

where the geopotential, U, is expressed in the Pines formulation using direction 
cosines, which is discussed in Chapter 2 and represented by Equation (2.4). The 
OSU86F gravity model (discussed in Chapter 2) through degree and order 5 was 
used for the geopotential coefficients in the above calculations. 

4.4.2 Drag 


The drag force per unit mass induced on the outer satellite is 



where C D is the outer satellite's coefficient of drag (C D = 3.5), S is frontal area of 
the spacecraft (S = 1.06 m 2 ), p is the density of the atmosphere, m is the mass of 
the satellite and V is the velocity of the spacecraft with respect to the atmosphere. 
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r 

pm 


coxr 


pm 


(4.9) 


The density of the Earth's atmosphere is much less predictable at altitudes 
above 125 km than at lower altitudes. The behavior of the atmosphere above 125 
km is due to many complex processes, such as geomagnetic heating, sunspot 
activity (solar flux), seasonal latitudinal variation of helium and the lower 
troposphere, and gravitational attraction of Earth's geoid. These processes 
determine the densities of the atmospheric constituents nitrogen (N 2 ), argon (Ar), 
helium (He), oxygen (O 2 and O), and hydrogen (H) which are involved in 
producing the exospheric temperature. 

The atmospheric density model used in this study was the Jacchia 1971 
Atmosphere Model with modifications by Roberts [ Jacchia , 1971]. This model 
calculates the densties of the atmosphere based on the aforementioned processes 
that raise or lower the exospheric temperature. Also included in this model are the 
analytical functions of the diurnal variation of the atmospheric bulge, the solar 
radiation flux at the 10.7 cm wavelength, the geomagnetic flux and the seasonal 
variations of the atmosphere. Since solar sunspot activity is known to occur with a 
1 1 year cycle, the 10.7 cm solar flux and the geomagnetic flux during the actual 
GRM mission (planned for January - June, 1990's) are predicted by using the data 



obtained from 1975 to 1986. The sunspot cycle's influence on the 10.7 cm solar 
flux and the geomagnetic planetary index for this period are presented in Figures 
4.6 and 4.7, respectively. 

The maximum sunspot activity occurred from 1980 to 1982, while the 
minimum sunspot activity occunred from 1985 to 1986. Atmospheric densities were 
calculated during these times to estimate the maximum and minimum densities the 
spacecraft would encounter during its mission in the 1990’s. The predicted 
latitudinal atmospheric density variation for altitudes of 150 to 170 km during the 
vernal equinox of 1991 (high solar activity) and the summer solistice of 1996 (low 
solar activity) are shown in Figures 4.8 and 4.9, respectively. Here it is found that 
the density of the atmosphere at 160 km altitude varies between 6 to 17 x 10* 10 
kg/m 3 for high solar activity and between 3.5 to 11.5 x 10* 10 kg/m 3 during the. 
minimum activity. 

Radial, along-track and cross-track drag components along the trajectory of 
the proof mass's 160 km high, drag-free orbit were computed using the density 
ranges stated previously for high and low solar activity; and they would resemble 
those shown in Figures 4.10 through 4.15. Since, at satellite altitudes, the drag 
induced on the spacecraft is porportional to the inverse of its mass, the drag will 
slowly increase during the lifetime of the mission as the fuel is expended. The drag 
values during the beginning (when the fuel tanks are full), the middle (half-full), 
and the end of the mission (empty) are also shown on the aforementioned figures. 

Along-track drag per unit satellite mass range from -30 to -180 microns/sec 2 
during high solar activity and from -20 to -140 microns/sec 2 during low solar 
activity. The unique patterns in these plots arise from the location on Earth that the 
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orbit passes. The polar orbit begins at the equator with a northemly direction; the 
peaks in the radial and along-track figures indicate the north and south pole regions 
where the magnitudes of drag are a minimum; the depressions indicate the 
equatorial regions where the magnitude of drag is greatest. The cross-track drag in 
Figures 4.14 & 4.15 arise from the rotation of the Earth and its atmosphere; the 
spacecraft's ascending track induces drag on its port side while the descending track 
induces drag on its' starboard side; and cross-track drag is zero at the poles. 

4.4.3 Thrust 


The thrust per unit mass generated by the thrusters is modelled as an on-off 
step function. When fore or aft thrusters need to fire then. 


thrust 


+16 Nt 

aft thrusters on 


on7 

off 

(4.10) 

-16N7 

fore thrusters on 



If cross-track thrusters need to fire then 


f 


thrust 


1 + 4 N n starboard thrusters on 
ONn off 

- 4 N n port thrusters on 


(4.11) 


or if radial thrusters need to fire then, 
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f =< 

thrust 


| + 4Nr top thrusters on 

0 Nr off (4.12) 

| - 4Nr bottom thrusters on 
where r, I, n denotes the radial, tangential (along-track) and normal (cross-track) 
directions, respectively, in a body-fixed RTN coordinate system. 


4.5 Numerical Integration of the Equations of Motion 


The Runga-Kutta Fehlberg 4(5) numerical integrator [Bettis, 1977] was used 
to integrate Equation 4.6 to determine the ephemeris of the proof mass's state vector 
with respect to the outer satellite and Equation 4.4 to determine the ephemeris of the 
proof mass. This fifth order, singlestep integrator allows for the integration of 
discontinuous functions such as the on-off switching of the thrusters. Embedded 
within this integrator is a fourth order method for automatic stepsize selection, but it 
is not used in this study. 

The step size (h) was set at 0.5 seconds during coasting and 20 milliseconds 
during thrusting phases of the proof mass trajectory with respect to the outer mass. 
The step size is further reduced in both phases to reach switch boundaries precisely. 

4,6 The Logic of the Thruster Control Algorithm 

The thruster control algorithm was programmed in FORTRAN-77 using UT 
System CHPC's Vax 8600 front-end computer and compiled and executed on the 
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Cray X-MP/24 Supercomputer. Conditional IF statements simulate the control 
system by monitoring the proof mass's position relative to the outer satellite in an 
RTN coordinate system and determining when, for how long, and what thrusters 
need to fire. 

The thruster control algorithm is implemented with the initial conditions and 
the various parameters listed in Tables 4.1 through 4.4. Starting at the beginning of 
the coast phase of the limit cycle with the numerical integration step size of 0.5 sec, 
the equations of motion of the proof mass with respect to the inertial frame and with 
respect to the centroid of the spacecraft (Equations 4.4 and 4.6) are integrated. At 
the end of every third step, the along-track position of the proofmass with respect to 
the centroid of the spacecraft and time is recorded. Upon recording the third data 
point, a parabola is fitted through the along-track position versus time curve and the 
along-track velocities of the ball with respect to the spacecraft's centroid is then 
calculated at these three positions. A parabola is fitted through the velocity versus 
position curve. Then the time of the aft thruster turn-on switch (deadband region) 
can be predicted at the intersection of this coast arc curve and an analytical 
representation of the velocity versus position curve (Equation 4.17) for the thrust 
phase of the limit cycle. If the difference of this time and the current time is less 
than the integration step size, the step size is reduced to this difference. 

Once the aft thruster tum-on switch has been activated, the drag force per unit 
mass is estimated using the total time of coast and the thruster-on and -off switch 
velocities (Equation 4.20). The aft thrusters are turned on and the integration step 
size is changed to 20 milliseconds. While integrating the ball’s state with respect to 
the spacecraft's centroid, the time of the thruster turn-off switch is predicted using 
the tum-on and turn-off velocities (explained below) and the estimated drag and 
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thrust per unit mass (Equation 4.21). Then, if the difference between the current 
time and the predicted time of thruster turn-off is less then the integration step size, 
the step size is reduced to this difference. Upon reaching the turn-off time, the aft 
thrusters are turned off and the limit cycle is repeated. 

Since the drag and thrust per unit mass changes according to solar activity, 
location above Earth and fuel expenditure, the imparted velocity at the thruster turn- 
off switch must be changed throughout the mission in order to keep the ball 
excursions within the designed -0.5 mm limit. When this value is held constant 
throughout the mission, either too small or too large limit cycles result which cause 
either the aft thrusters to fire very frequently or the fore thrusters to fire, thereby 
greatly increasing fuel consumption. Therefore, this imparted velocity turn-off 
switch was determined for each phase (beginning, middle and end) of the mission 
such that the maximum ball along-track excursion remains within the designed -0.5 
mm limit, but was held constant in each case. These values are listed in Table 4.4 
for high and low solar activity, respectively, and they range from -0.44 to -0.24 
mm/sec. 

Unless unexpected disturbances cause the proof mass to reach the radial or 
cross- track ±1.0 mm deadband regions, there is no need to fire radial or cross- 
track thrusters. Therefore, the computational thruster control algorithm outlined 
below involves only along-track drag compensation. All radial and cross-track 
excursions are controlled by the ball-centering attitude control laws explained in 
Section 4.7. 
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Along-track Drag-free Thruster Control Algorithm 

1. Specify: a) Initial Conditions of proof mass state with respect to the outer 
mass and the proof mass state in the inertial mean of 2000 coordinate frame 

b) Thrust model parameters, 

c) Gravity model parameters, 

d) Drag model parameters, and 

e) Control law parameters 

The corresponding parameters used in this study are listed in Table 4.1 through 4.4 
and Equation 4.10. 

2. Set step-size to 0.5 seconds. 

3. Integrate equations of motion Equations 4.4 and 4.6. 

4. Record time and proof mass position with respect to outer mass every third 

step. 

5. Once three data points of time and along-track position have been recorded 
along the coast arc phase, the along-track control law is implemented by fitting a 
parabola of the form 

x = t 2 a : + t a 2 + 83 (4.14) 

through the curve using a least squares approximation where a lf a 2 , a 3 are the 
coeffients that describe the curve, x is the along-track position and t is the time since 
the coasting phase began. 

6. Determine velocities at these three data points, 

x. = 2t i a 1 + a 2 i = 1 , 2,3 (4.15) 

7. A parabola is fitted through the three position and velocity data points to 
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find the coefficients and b 3 of the equation 

x = x 2 ^ + xb 2 + b 3 (4.16) 

This equation is used to describe the phase plane characteristics of the coast arc in 
the along-track direction. 

8. The analytical representation of the phase plane motion in the along-track 
direction during the thrust arc is, 

I i2 = -< f D,» g + f Tta» )(X - (4.17) 


where x^ is the maximum excursion of the proof mass with respect to the outer 
satellite (taken to be +1.0 mm), frhrust the thrust per unit satellite mass fired in the 
transverse direction, and f^g is the estimated transverse drag force per unit mass 
(initially set at a nominal value to initiate first control law predictions). 

9. The velocity at the transition point from the coast arc to the thrust arc is 
defined by the intersection of Equations 4.16 and 4.17 and is predicted using the 
quadractic equation (positive root), 


-B ±J B 2 - 4AC 
*ON “ — 2A 


(4.18) 


where, 


A = bj 
B = b 2 
C = hj 


1 

^Thrust + ^Drag) 


X 


max 


(4.19) 


The position and time of deadband intersection (thruster tum-on switch, Xq N , tQ N ) 
is calculated using Equations 4.16 and 4.15 respectively. 
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10. If the difference between the predicted turn-on switch time tQ N and the 
current time is less than the integration step size, the step size is reduced to this 
difference. The time of entire coast phase, tQ^, is recorded. 

1 1. Estimate drag accelerations during the coast phase 


Drag 


X OFF ■ X ON 


Coast 


(4.20) 


where XQppis the velocity at the thruster turn-off switch (thrust imparted velocity of 
the proof mass relative to the spacecraft) calculated using Equation 4.17 where x^ 
is the maximum excursion of the ball when the magnitude of drag is a minimum. 

12. Tum-on transverse thrusters, change integration step size to 20 
milliseconds. 

13. Thruster turn-off time, toFF is then predicted using 


_ x on ' X OFF 
l OFF _(£ 


- f ) 

Thrust Drag'' 


(4.21) 


14. If the difference between the current time along the thrust arc and toFF is 
less than the integration step size, reduce the step size to this difference to reach the 
turn-off switch precisely. 

15. Once t = toFF> the thrusters are turned off and steps 2 - 15 are repeated. 
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Table 4.1 Initial Conditions 

The initial conditions of the proof mass's state in inertial geocentric equatorial mean of 
1950 system are: 


X = 

6538137.0 m, 

v x = 

0.001 m/s, 

Y = 

0.001 m, 

V Y = 

0.001 m/s, 

Z = 

0.001 m, 

V Z = 

7808.03729 ml 

the initial conditions of the proof mass's state with respect to the outer satellite are: 

e R = 

0.00 mm, 

e R = 

0.00 mm/s, 

e T = 

0.95 mm, 

e T = 

Xqff mm ^ s 




(see Table 4.4) 

e N = 

0.00 mm, 

£ N = 

0.00 mm/s. 

Table 4.2 Parameters for the Gravity Model 


Earth's mean radius, 


Ae = 6378137.0 

m, 


Earth's gravitational parameter, ft = 3.9860044 x 10 + * 4 m-Vs^, 

Earth's rotational velocity, to = 7.29211611 x 10*® rad/s, 

Gravity Coefficients, Cnm & Snm from 5x5 OSU86F gravity field 
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Table 4.3 Parameters for the Drag Model 


Mass of satellite: 

m 

= 

1334 kg, 

Mass of Fuel: 


= 

1400 kg 

Total Mass at launch from space shuttle, 


= 

2734 kg 

at mission start (95% Full), 


= 

2664 kg, 

middle of mission (50% Full), 


= 

2034 kg 

near end of mission (5% Full), 


= 

1404 kg, 

Satellite's projected cross-sectional area. 

S 

= 

1.06 m 2 , 

Satellite's coefficient of drag, 

c D 

= 

3.5, 


Tahle 4.4 Parameters for the Control Laws 


The maximum forward designed excursion, X niax = 

+1.0 mm, 



The minimum aft designed excursion, X m j n = 

-0.5 mm, 



• 

The thruster off switch, Xopp, 




95% Full 

S0% Full 

5% Full 

High Solar Activity: -0.308 mm/s 

-0.352 mm/s 

■0.437 

mm/s 

Low Solar Activity: -0.240 mm/s 

-0.274 mm/s 

-0.330 

mm/s 


4.7 Ball-Centering Attitude Control Law 

A ball-centering attitude control law is coupled with the DISCOS translational 
control system so that radial and cross-track excursions of the proof mass with 
respect to the outer mass can be stopped without firing the radial and cross-track 
thrusters. Figure 4.16 shows how a small pitch bias, a, about the cross-track axis. 


produces a small aerodynamic lift force. 


*0^ - f Dr a&r a 


(4.21) 


Likewise, a small aerodynamic sideslip force can be produced to center the proof 
mass in the cross-track direction by introducing a small yaw bias, (3, about the 
radial axis (Figure 4.17), 


f Sideslip N - W&jP (4.22) 

These angles are computed by the ball-centering attitude feedback control laws of 
the form 

a = - K i n x r ' ^ x r (4.23) 

p = K 1 r*N + (4.24) 

where K 1r and K 2r and K 1n and K 2n are the feedback gains about the radial and 
normal axis respectively (listed in Table 4.5). The gains, K 1n and K 2]s} , were 
determined by Ray and Jenkins [1982], and the K 1r and K 2r gains were 
determined in this study through an iterative process. 


Table 4.5 — Ball-Centering Attitude Control Gains 
Pitch gains (about the normal axis), K 1n = 3.0, 

K 2n = 200.0, 

Yaw gains (about the radial axis), k 1r = 100 -°’ 

K *r * 


3000.0. 
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4.8 Results 

The GRM DISCOS simulation using the drag-free thruster control algorithm 
was performed during high and low solar activity and during the beginning (95% 
full), middle (50% full), and end (5% full) of mission to estimate fuel expenditures 
during the mission. Each simulation of one orbital revolution required 
approximately 86 seconds of CPU time on the Cray X-MP/24 supercomputer. The 
number of thruster firings and fuel expended per orbit and 6 month mission along 
with other pertinent data are presented in Table 4.6 and 4.7 respectively for the high 
and low solar activity simulations. 

■4.8.1. The Effectiveness of the Drag-Free Thruster Control Algorithm 

Along-track proof mass excursions with respect to the centroid of the outer 
mass against time for one orbit during high and low solar activity are shown in 
Figures 4.18-a, -b, and -c and Figures 4.19-a, -b, and -c when the spacecraft has 
respectively, 95% fuel, 50% fuel and 5% fuel of initial fuel onboard. Radial and 
cross-track proof mass excursions are shown in Figures 4.20-a, -b, and -c for high 
solar activity and Figures 4.21-a, -b, and -c for low solar activity. These excursions 
were plotted during two orbits to verify that they were not securly increasing. 
Along-track phase plane diagrams in Figures 4.22-a, -b, and -c and Figures 4.23-a, 
-b, and -c show the limit cycles that are produced from the control laws respectively 
during the high and low solar activities. 

The aforementioned figures demonstrate the ability of the thruster control 
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algorithm to offset effects due to drag and keep the GRM spacecraft essentially 
drag-free. The along-track phase plane diagrams show the clock-wise pattern of the 
limit control cycles where the nearly straight (actually parabolic) portion of the cycle 
at the +1.0 mm region is the thrusting arc from top to bottom and the band of 
parabolas with their apexes of +0.5 mm to -0.5 mm are the coast arcs from bottom 
to top. They show that when the proof mass reaches the aft deadband limit of + 1.0 
mm inside the spacecraft’s cavity with a positive approaching velocity, the along- 
track thrusters on the spacecraft fire from 55 to 105 milliseconds until a 
predetermined negative imparted velocity is achieved. Then the proof mass is seen 
coasting until a maximum excursion of -0.5mm relative to the spacecraft's cavity is 
reached at zero velocity. Drag then overcomes the thrusting effort and causes the 
proof mass to "fall" with positive velocity back to the deadband region. 

4.8.2 The Effectiveness of the Ball- Centering Attitude Control Law 

Figures 4.20 and 4.21 (-a, -b, and -c) and 4.24 to 4.28 demonstrate the 
effectiveness of the ball-centering attitude control law to produce small lift and 
sideslip forces to center the proof mass in the radial and cross-track directions. 
Figures 4.20(-a, -b, and -c) and 4.21(-a, -b, and -c) show that the radial proof 
mass excursions deviate no more than +0.40 mm or no less than -0.32 mm from 
the centroid of the spacecraft. These figures also show that the cross- track ball 
excursions deviate between ±0.60 mm.. Radial and cross-track relative velocities as 
a function of time are plotted along two orbits in Figure 4.24 during high solar 
activity and 50% fuel and are indicative of all high and low solar activity cases. The 
radial velocities range from -1.5 to +1.5 jim/s while the cross-track velocities range 
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from -0.75 to +0.75 (im/s. Also indicative of all cases, the ball-centering attitude 
control pitch and yaw angles for the above case are shown in Figures 4.25 and 4.26 
during two orbits. It is shown that very small pitch angles of 0.057 to 0.069 
degrees are produced to counteract radial excursions while relatively large yaw 
angles of ±3.4 degrees are needed to counteract the cross-track excursions. Finally, 
radial and cross-track phase plane diagrams for this case are shown in Figures 4.27 
and 4.28. The irregularities at the beginning of each plot, especially in the cross- 
track curves, resulted from the initial guess of the magnitude of drag to initiate the 
algorithm. Once initiated the curves become periodic with the orbit; local minimum 
and maximum excursions are indicative of the drag induced on the spacecraft in 
these directions or its location above the Earth when compared to the drag profiles 
in Figures 4.10 through 4.15. 

Without the ball-centering attitude control law, the radial and cross-track 
excursions would hit the ±1.0 mm deadband regions causing both radial and cross- 
track thrusters to fire briefly every 100 to 200 seconds. 

4.8.3 Imparted Velocity Thurster Tum-off Adjustments 

The along-track phase plane diagrams in Figures 4.22 and 4.23 (-a, -b, -c) 
show the results of the velocity tum-off switch adjustments. Each diagram shows 
the maximum excursion to be approximately -0.5 mm, resulting in maximum coast 
arc times of 14 to 25 seconds (as indicated in Tables 4.6 and 4.7) depending on 
solar activity and the amount of fuel onboard. However, holding the imparted 
velocity tum-off switch constant through the drag variation of one orbit results in 
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m i nim um excursions of +0.50 mm for high solar activity and +0.60 mm for low 
activity. These minimum excursions result in relatively short coast times between 
5.0 and 7.5 seconds (as indicated in Tables 4.6 and 4.7) which may not be 
acceptable for minimizing fuel consumption or for post ball reconstruction. 
Therefore, the control algorithm should be able to adaptively adjust the tum-off 
switch due to these effects. However, the selection of an adaptive control law is 
beyond the scope of this research. 

4.8.4 Fuel Expenditure 

Both drag and thrust per unit mass increase when fuel is expended because of 
the decrease in spacecraft mass. For this reason, three senarios during the mission 
are studied: 1) the beginning of the mission (after orbit injection) when the 
spacecraft has approximately 95% fuel left onboard, 2) the middle of the mission 
with 50% fuel onboard, and 3) near the end of the mission when there is 
approximately 5% fuel left onboard. Since drag increases during high solar activity, 
fuel expenditure is expected to be much greater during this period than during low 
solar activity. To best estimate fuel expenditures, both high and low cases must be 
investigated. 

The fuel expenditure for the GRM mission is estimated by first determining 
the mass fuel rate of the thrusters and the total thruster-on time during one complete 
orbital revolution about Earth. Then assuming that fuel consumption will be nearly 
constant through the entire six month mission, the fuel expended per orbit is 
multiplied by the number of orbital revolutions at the end of six months. 

The mass fuel rate of the thrusters when fired is 
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m c . = 

Fuel 


' Thrust 


Se 1 ; 


SP 


(4.25) 


where F-pj^st is the total thrust produced by the aft thrusters, g e is the Earth's 
gravity at the surface and Isp is the specific impulse of the fuel. With an I SP of 200 
seconds and F ThnJst of 16 Newtons (g e = 9.81 m/s 2 ), the fuel mass rate is 


0.0081633 kg/s. 

The fuel expenditure for the entire six month mission can be estimated by 
using the above algorithm to find the total thruster-on time, k)N Tota i’ f° r one orbital 
revolution about Earth. The fuel expended per orbit, is then 


Fuel ^ 

ObiT " m Fud t ONrotal 


(4.26) 


and the mission fuel expenditure can be estimated by noting that the six month 
mission requires approximately 2953 revolutions about Earth. 

As listed in Tables 4.6 and 4.7, thruster on-times varied between 75.5 to 105 
milliseconds during high solar activity, and between 58.9 to 81.8 milliseconds 
during low activity. The number of firings per orbit rose from 423 during the 
beginning of the mission to 583 firings near the end for the high solar activity case 
and from 395 to 543 firings during the low activity. However, the total on-time per 
orbit from beginning to end of the mission was nearly constant at approximately 
43.8 sec during high solar activity and 31.7 sec for low activity. The idea that the 
fuel expended per orbit is independent of the amount of fuel left in the tanks was 
proven by observing the values of fuel expended per orbit listed in Tables 4.6 and 
4.7. Here, the fuel expended per orbit was also nearly constant at 0.357 kg/orbit 
during high solar activity and 0.256 kg/orbit for low activity. Finally, with 
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approximately 2953 orbital revolutions, fuel expenditure estimates of 763 and 1056 
kilograms of hydrazine (as listed in Tables 4.6 and 4.7) were determined, 
respectively during the low and high solar activities. 

These fuel expenditures were compared with other published GRM fuel 
expenditure estimates. According to Ray and Jenkins [1982], their DISCOS 
simulation resulted in a cycle time of 10-20 seconds and about two million firings 
for the six month mission. Using a totally diffuse reflection drag model for worst 
case estimates, their simulations estimated the fuel expenditure for a seven month 
mission to be approximately 1150 kg giving 2.6 million N-sec of total impulse to 
counteract drag. This corresponds to nearly 986 kg of hydrazine for a six month 
mission if the fuel consumption rate was constant during entire mission. Their 
optimistic fuel estimate used only purely specular reflection and required 660 kg for 
seven months which would correspond to 566 kg for six months. According to 
Keating et al, [1986], the Goddard Space Flight Center's Flight Dynamics 
Division determined the 180 day (six month) mission would require approximately 
1051 kg of hydrazine using worst case assumptions such as high solar activity. 

Table 4.8 lists the worst and best case fuel expenditures of each study. It 
shows that the results of this report were reasonably close when compared to the 
other studies; the results were 0.5% greater than Keating et al, [1986] and 7.0% 
greater than Ray and Jenkins [1982] for the worst case, and 35% greater than Ray 
and Jenkins [1982] for the best case. 
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Results of the GRM DISCOS Simulation during High 
95% Full 126184 kg) M (2034 kg) 


IMl AA 


Imparted velocity 

No. of firings/orbit 

No. of firings/6 month mission 

Fuel expended/orbit 

Six month fuel expenditure 

Maximum thrust time 

Minimum thrust time 

Maximum coast time 

Minimum coast time 

Total thrust time/orbit 

Maximum estimated drag magnitude 

Minimum estimated drag magnitude 


■0.30795 mm/s 

•0.35244 mm/s 

423 

484 

1.25 million 

1.43 million 

0.35749 kg 

0.35741 kg 

1055.7 kg 

1055.5 kg 

105.01 msec 

91.76 msec 

104.06 msec 

90.93 msec 

19.392 sec 

16.942 sec 

7.591 sec 

6.248 sec 

43.793 sec 

43.783 sec 

86.140 pm/s 2 

112.82 pm/s 2 

31.781 pm/s 2 

41.624 pm/s 2 


Solar Activity 
5%_Full (1404 kel 
-0.43695 mm/s 
583 

1.72 million 
0.35767 kg 
1056.3 kg 
76.24 msec 
75.53 msec 
14.075 sec 
5.347 sec 
43.815 sec 
163.44 pm/s 2 
60.302 pm/s 2 
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Table 4.7 Results of the GRM DISCOS Simulation during Low Solar Activity 


95% Full (2664 kg) 

50% Full (2034 kg) 

5% Full (1404 kgi 

Imparted velocity 

•0.23957 mm/s 

■0.27417 mm/s 

-0.32999 mm/s 

No. of firings/orbit 

395 

451 

543 

No. of firings/6 month mission 

1.17 million 

1.33 million 

1.60 million 

Fuel expended/orbit 

0.25839 kg 

0.25850 kg 

0.25857 kg 

Six month fuel expenditure 

763.07 kg 

763.39 kg 

763.60 kg 

Maximum thrust time 

81.84 msec 

71.51 msec 

59.42 msec 

Minimum thrust time 

81.19 msec 

69.94 msec 

58.94 msec 

Maximum coast time 

24.959 sec 

21.813 sec 

18.134 sec 

Minimum coast time 

7.058 sec 

6.167 sec 

5.123 sec 

Total thrust time/orbit 

31.653 sec 

31.667 sec 

31.675 sec 

Maximum estimated drag magnitude 

67.889 pm/s 2 

88.916 pm/s 2 

128.81 pm/s 2 

Minimum estimated drag magnitude 

19.192 pm/s 2 

25.136 pm/s 2 

36.416 pm/s 2 


Table .4,8 Comparisons of Fuel Expenditure s for 6 Mon th Mission 


Worst Case 


Keating et al. T19861 Ray & Jenkings. T19811 
1051 kg 986 kg 


N.A. 


566 kg 


Best Case 


This Report 
1056 kg 
763 kg 


PRE-AMP ELECTRONICS 
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Figure 4.1 The GRM Disturbance Compensation System (DISCOS) Mechanism. 
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SST One-Way Doppler Antennae (42 GHz) 
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Figure 4.5 The Vector representing the Proof Mass 
with respect to the Outer Satellite 


SOLAR FLUX 


Figure 4.6 Solar flux data in 10e-22 
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Figure 4.10 Predicted drag profiles along 2 GRM (160km) orbits on the 
vernal equinox of 1991 (high solar activity). 
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4.1 1 Predicted drag profiles along 2 GRM (160km) orbits on the 
summer solistice of 1996 (low solar activity). 
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Figure 4.12 Predicted drag profiles along 2 GRM (160km) orbits on the 


vernal equinox of 1991 (high solar activity). 
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Figure 4.13 Predicted drag profiles along 2 GRM (160km) orbits on the 


summer solistice of 1996 (low solar activity). 
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Figure 4.14 Predicted drag profiles along 2 GRM (160km) orbits on the 
vernal equinox of 1991 (high solar activity). 
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Figure 4.15 Predicted drag profiles along 2 GRM (160km) orbits on the 
summer solistice of 1996 (low solar activity). 
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Figure 4.17 Aerodynamic sideslip force generated by 
the Ball Centering Control Law. 
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Figure 4.18 : a Along-track position of proof mass with respect to spacecraft’s 
centroid against time for one orbit during high solar activity with 95% fuel onboard. 
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Figure 4.18-b Along-track position of proof mass with respect to spacecraft's 
centroid against time for one orbit during high solar activity with 50% fuel onboard. 
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Figure 4.18-c Along-track position of proof mass with respect to spacecraft's 
centroid against time for one orbit during high solar activity with 5% fuel onboard. 
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Eigprg 4, 19- a Along-track position of proof mass with respect to spacecraft's 
centroid against time for one orbit during low solar activity with 95% fuel onboard. 
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Figure 4.19-b Along-track position of proof mass with respect to spacecraft's 
centroid against time for one orbit during low solar activity with 50% fuel onboard. 
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Figure 4.19-c Along-track position of proof mass with respect to spacecraft's 
centroid against time for one orbit during low solar activity with 5% fuel onboard. 
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Figure 4.20-a Radial and cross-track position of proof mass with respect to 
spacecraft's centroid against time for two orbits during high solar activity with 95% 

fuel onboard. 
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Figure 4.20-b Radial and cross-track position of proof mass with respect to 
spacecraft's centroid against time for two orbits during high solar activity with 50% 


fuel onboard. 
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Figure 4.20-c Radial and cross-track position of proof mass with respect to 
spacecraft's centroid against time for two orbits during high solar activity with 5% 

fuel onboard. 
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Figure 4.21 -a Radial and cross-track position of proof mass with respect to 
spacecraft's centroid against time for two orbits during low solar activity with 95% 

fuel onboard. 


96 


O 

« 



TIME (MIN) 


Figure 4.21-b Radial and cross- track position of proof mass with respect to 
spacecraft's centroid against time for two orbits during low solar activity with 50% 

fuel onboard. 
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Figure 4.2 1-c Radial and cross-track position of proof mass with respect to 
spacecraft's centroid against time for two orbits during low solar activity with 5% 

fuel onboard. 
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Along-track phase plane limit cycle for one orbit during high solar 
activity with 50% fuel onboard. 
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Figure 4.22-c Along-track phase plane limit cycle for one orbit during high solar 

activity with 5% fuel onboard. 
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Figure 4.23-c Along-track phase plane limit cycle for one orbit during low solar 

activity with 5% fuel onboard. 
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Figure 4.24 Radial and cross-track proof mass velocities with respect to the 
spacecraft against time for two orbits during high solar activity with 50% fuel 

onboard. 
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Figure 4.25 Ball-centering attitude control pitch angle against time for two orbits 
during high solar activity with 50% fuel onboard. 
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Figure 4.27 Radial phase plane diagram (velocity vs. position) for one orbit dining 
high solar activity with 95% fuel onboard. 
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CHAPTER V 

CONCLUSIONS 

Presented in the first part of this study are two orbit simulations, one 
representing the actual GRM orbits and the other representing the orbit estimated 
from orbit determination techniques. These orbit simulations are essential for the 
geodetic science community to test gravity evaluation techniques before the GRM 
mission is launched. The second part of this study involved creating a computer 
algorithm to simulate GRM's drag compensation mechanism (DISCOS) so that fuel 
expenditure and proof mass trajectories relative to the spacecraft centroid could be 
calculated for the mission. 

Although the low-low GRM spacecraft configuration was used exclusively, 
the orbit information obtained in this study is applicable to all possible scenarios 
since at least one low satellite will be used. In addition, the drag compensation 
study is also still applicable to a possible change of the design of DISCOS 
mechanism in the low-low scenario to a "two- stage" design, or the replacement of 
the SST (and elimination of the second satellite) with a gravity gradiometer. 

The two-stage DISCOS and gravity gradiometer take the single stage DISCOS 
and SST a step further by eliminating the need for post flight proof mass velocity 
reconstruction since measurements taken by these instruments are always taken 
about the proof mass centroid. The two- stage DISCOS system incorporates an 
intermediate stage that is equipped with the SST antennae and is magnetically forced 
to follow the proof mass constantly while the outer stage uses the same concept as 
performed in this study to follow the intermediate stage by firing appropriate 
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thrusters to compensate for drag. The gravity gradiometer uses a set of high 
precision accelerometers whose measurements would replace the range rate 
measurements from the SST for the determination of the higher degree and order 
terms of the geopotential. 

5 t l_ S um m ary 

With a geopotential of degree and order 360, the 32 day true orbit simulation 
of the two GRM spacecrafts' proof masses required a total of 19.2 hours of CPU 
time on the Cray X-MP/24 supercomputer, nearly 100 times faster than on any 
conventional mainframes computers. The initial conditions determined by White 
[1987] satisfied the requirement that the orbit repeat after exactly 32 sidereal days 

The orbit determination solution discussed here merely indicates some of the 
major points that must be considered when solving for the nominal trajectory. It 
was shown that the orbit can be determined within the 3 a accuracy requirement of 
100 meters in the radial and 300 meters in the along- and cross-track directions for 
the gravity mission. The magnetic mission requirements of 60 meters radial and 100 
meters in the along- and cross-track components could not be met with this 
solution; however the nominal orbit simulation will be used to further improve 
gravity models so that eventually the residuals will meet the magnetic mission 
requirements. 

The one-way Doppler algorithm, created for the SST simulation, rarely did 
not converge for measurements received by satellite two's antenna, and the latest 
measurements calculated during these non-convergences did not meet the 1 (im/s 
requirement. Much effort went into finding the cause of this iteration problem, but 
the problem was never solved completely. It was finally thought to be due to 
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round-off errors occurring during the interpolation of the ephemeris between the 
recorded time points. 

The results of the GRM DISCOS simulation demonstrated that the spacecraft 
can essentially be "drag-free", that is, the spacecraft can be made to follow the 
purely gravitational trajectory of the proof mass. The results showed that the 
centroid of the spacecraft can be controlled so that it will not deviate more than 1.0 
mm in any direction from the centroid of the proof mass. 

Radial and cross-track proof mass excursions were shown to be constrained 
within the ±1.0 mm region without firing radial or cross-track thrusters, but with 
the use of the ball-centering attitude control law which pitches or yaws the 
spacecraft at slight angles to produce counteracting lift or sideslip forces. 

The aft along-track thrusters were determined to fire for 55 to 105 

* 

milliseconds every 5 to 25 seconds depending on drag magnitudes and the current 
satellite mass. The total number of aft thruster firings ranged from 1.2 to 1.7 
million. Total thrust-on times, and thus fuel expenditure rates, were essentially 
independent of the variation in mass due to fuel consumption for either the high 
solar or low solar activity case. It was also shown that estimated fuel expenditures 
of 763 to 1056 kg of hydrazine for drag compensation compared closely to other 
studies. 

5.2 Future Study 

Now that both a representation of the proof mass orbit trajectories and a way 
of simulating the drag compensation mechanism exists, the simulation of the outer 
satellites' orbits can be determined. Proof mass displacements from the satellites' 
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mass centers corrupt the SST range-rate measurements, thus the outer satellite orbit 
simulation will incorporate this source of error thereby allowing precise calculations 
of the SST range-rate measurements. 

Many refinements and improvements are needed to improve the drag-free 
thruster control algorithm. To fully simulate the DISCOS mechanism, adaptive 
control laws should be implemented to provide coasting arcs of nearly constant 
duration. An algorithm to adjust the thruster turn-off switches (imparted relative 
velocities) according to the magnitude of drag and satellite mass is needed since 
either too small or too large limit cycles would result if this parameter is held 
constant throughout the mission. This would require aft thrusters to fire frequently 
or fore thrusters to fire occasionally, thus increasing fuel consumption. It was not 
clear from previous reports that this point has been studied. 

The Jacchia 1971 atmospheric model used in this simulation was also slightly 
outdated. A more current one, such as the Jacchia 1977 atmospheric model should 
be incorporated into the study. In addition, the drag model used in this study did 
not take into account the effects due to complicated atmospheric gas impingements 
and reflections. Also, a detailed model of the satellite’s surface area incident to the 
airstream is needed since the incidental area (frontal area) used in this study was 
held constant and radial and cross-track drag components were determined from 
this value. 

Finally, the inclusion of a fuel expenditure algorithm into the DISCOS 
simulation is needed to better calculate fuel expenditures. Longer simulations 
(greater than one orbit ) would further improve these calculations. 
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